{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A time-dependent problem, Burgers' equation\n",
    "\n",
    "We will solve the viscous Burgers equation, a nonlinear equation for the advection and diffusion on momentum in one dimension.\n",
    "\n",
    "$$\n",
    "\\frac{\\partial u}{\\partial t} + u \\frac{\\partial u}{\\partial x} - \\nu \\frac{\\partial^2 u}{\\partial x^2} = 0\n",
    "$$\n",
    "\n",
    "We will solve on a periodic interval mesh, and therefore do not impose any boundary conditions.  As usual, we need to derive a variational form.\n",
    "\n",
    "## Spatial discretisation\n",
    "\n",
    "We first discretise in space, mulitplying by a test function $v \\in V$ and integrating the viscosity term by parts to obtain the semi-discrete problem. Find $u(x, t) \\in V$ such that\n",
    "\n",
    "$$\n",
    "\\int_\\Omega \\frac{\\partial u}{\\partial t} v + u \\frac{\\partial u}{\\partial x} v + \\nu \\frac{\\partial u}{\\partial x}\\frac{\\partial v}{\\partial x} \\, \\mathrm{d}x = 0 \\quad \\forall v \\in V.\n",
    "$$\n",
    "\n",
    "## Time discretisation\n",
    "We now need to discretise in time.  For simplicity, and stability we'll use backward Euler, replacing all instances of $u$ with $u^{n+1}$ and the time derivative by $\\frac{u^{n+1} - u^n}{\\Delta t}$.  We end up with the discrete problem, find $u^{n+1} \\in V$ such that\n",
    "\n",
    "$$\n",
    "\\int_\\Omega \\frac{u^{n+1} - u^n}{\\Delta t} v + u^{n+1} \\frac{\\partial u^{n+1}}{\\partial x} v + \\nu \\frac{\\partial u^{n+1}}{\\partial x}\\frac{\\partial v}{\\partial x} \\, \\mathrm{d}x = 0 \\quad \\forall v \\in V.\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "\n",
    "To solve the problem in a concrete setting, we need two things.  A domain, and an initial condition for $u$.  For the former, we'll choose a periodic interval of length 2, for the latter, we'll start with $u = \\sin(2 \\pi x)$.\n",
    "\n",
    "In addition we need to choose the viscosity, which we will set to a small constant value $\\nu = 10^{-2}$\n",
    "\n",
    "As ever, we begin by importing Firedrake"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from firedrake import *\n",
    "\n",
    "n = 100\n",
    "mesh = PeriodicIntervalMesh(n, length=2)\n",
    "\n",
    "x = SpatialCoordinate(mesh)[0]\n",
    "\n",
    "u_init = sin(2*pi*x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "nu = Constant(1e-2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We choose degree 2 piecewise continuous Lagrange polynomials for our solution and test space:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "V = FunctionSpace(mesh, \"Lagrange\", 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need solution functions for $u^{n+1}$ and $u^n$, along with a test function $v$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "u_n1 = Function(V, name=\"u^{n+1}\")\n",
    "u_n = Function(V, name=\"u^{n}\")\n",
    "v = TestFunction(V)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We provide the initial condition for $u_n$, and choose a $\\Delta t$ such that the advective Courant number is around 1.  This is more restrictive than required for stability of the time integration, but gives us enough accuracy to see the temporal evolution of the system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "u_n.interpolate(u_init)\n",
    "dt = 1.0 / n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we're ready to define the variational form.  Since this problem is nonlinear, note that we do not have a trial function anywhere.  We just write down the residual, Firedrake will automatically compute the Jacobian by differentiating the residual inside the nonlinear solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "F = (((u_n1 - u_n)/dt) * v +\n",
    "     u_n1 * u_n1.dx(0) * v + \n",
    "     nu*u_n1.dx(0)*v.dx(0))*dx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For visualisation purposes, we will save a copy of the state $u_n$ at each timestep, we can plot and animate these in the notebook if the `ipywidgets` package is installed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If passed an existing Function object, the Function \n",
    "# constructor makes a copy.\n",
    "results = [Function(u_n)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we loop over the timesteps, solving the equation and advancing in time.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = 0.0\n",
    "t_end = 0.5\n",
    "while t <= t_end:\n",
    "    solve(F == 0, u_n1)\n",
    "    u_n.assign(u_n1)\n",
    "    results.append(Function(u_n))\n",
    "    t += dt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This interactive plot should provide a slider that controls which iteration is plotted.  If you do not see it, it may be because you have not installed the necessary [ipython widgets](http://ipywidgets.readthedocs.io/en/latest/index.html).\n",
    "\n",
    "In the activated virtualenv in the terminal, you will need to do:\n",
    "\n",
    "```shell\n",
    "pip install ipywidgets\n",
    "jupyter nbextension enable --py widgetsnbextension --sys-prefix\n",
    "```\n",
    "\n",
    "and then restart the notebook kernel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU9b34/9c7G4EIkkhACEvCDoICRvZFASkUhVrRorctVi22vbXr1632tuqt9+LtvXXp7a91vdXW1irKYgUxLCICASKLgLKjhEXABBACIdv798ecM4whIcssZybzfj4e88jMWd8z+ZzzPp/P5yyiqhhjjIlfCV4HYIwxxluWCIwxJs5ZIjDGmDhnicAYY+KcJQJjjIlzSV4H0Bht2rTR7Oxsr8MwxpiY8sEHH3yuqpnVh8dkIsjOzqagoMDrMIwxJqaIyKc1DbemIWOMiXOWCIwxJs5ZIjDGmDhnicAYY+KcJQJjjIlzIUkEIvKCiBwRkS21jBcReUpEdonIhyIyKGDcDBHZ6bxmhCIeY4wx9ReqGsGfgYkXGD8J6OG8ZgJ/BBCRDODXwBBgMPBrEUkPUUzGGGPqISTXEajqeyKSfYFJpgIvqe+e1/ki0lpE2gNXA3mqWgwgInn4EsrfQxFXtDt27BirV69m7dq1lJeXc+WVVzJ48GA6duzodWimCSgtLWXNmjWsXLmSEydO0Lt3b0aPHk23bt28Ds1EmUhdUJYFFAZ83u8Mq234eURkJr7aBJ07dw5PlBFSVFTEAw88wEsvvQTA2bNnAWjWrBnl5eVMnjyZxx9/3DZY0yilpaXMmjWLJ598kjNnznypfJ09e5ahQ4fyu9/9jmHDhnkcqYkWMdNZrKrPqGququZmZp53hXTMePfddxkyZAjPPvssZWVldO3alV/84hf8+te/ZsCAASQlJfHmm28yaNAgXnnlFa/DNTHm448/Zvz48Tz88MMcP36czMxMfvKTn/Dv//7vjBw5klatWpGfn8+IESN47LHHqKqq8jpkEw1UNSQvIBvYUsu4p4FbAj5vB9oDtwBP1zZdba8rr7xSY9Ff//pXTUpKUkCHDRum27ZtO2+aw4cP6ze/+U0FFNBf/vKXWlVV5UG0JtasXLlS09PTFdD27dvr8uXLz5umpKREf/GLX/jL17Rp07S0tNSDaI0XgAKtaR9d08DGvOpIBJOBhYAAQ4G1zvAMYC+Q7rz2Ahl1rSsWE8GcOXNURBTQe+65RysrKy84/bPPPquJiYkK6H/8x39EKEoTqzZt2qStWrVSQKdMmaIlJSUXnH7p0qV68cUXK6DTp0+vszyapiGsiQBf5+4hoBxfO/8dwPeA7znjBfgDsBvYDOQGzHs7sMt5fac+64u1RLBp0yZNS0tTQB9++OF6z/faa6/5k8ecOXPCGKGJZYcPH9acnBwF9KabbtLy8vJ6zbd+/Xpt2bKlAvrII4+EOUoTDcJeI4jkK5YSwfHjxzU7O1sB/da3vtXgZp7f/e53CmhaWppu3749TFGaWFVRUaFjx45VQK+88ko9c+ZMg+Z/5513FNCEhAR96623whSliRaWCDwyffp0/0Z69uzZBs9fVVWlt956qwI6cODARi3DNF3/+Z//qYC2bdtWDx061KhlPPbYYwpoZmamHjx4MMQRmmhiicADc+fO9R/N79ixo9HLOXHihL/q/+ijj4YwQhPLtm3bpikpKQrowoULG72cyspKHT9+vAJ68803hzBCE20sEUTYyZMntVOnTgrok08+GfTyFi9erICmpqbqrl27QhChiWVVVVX+JqHbb7896OV98skn2qJFi6CTiolutSWCmLmOINb893//N4WFhQwaNIh//dd/DXp548aN41/+5V8oLS3lvvvuC0GEJpa9/vrrLFu2jEsuuYTHHnss6OV16dKFhx56CICf//znlJWVBb1ME0Nqyg7R/or2GkFhYaGmpqYqoCtWrAjZcg8ePOhfbk3niJv4cPbsWe3atasC+r//+78hW25ZWZm/CfKpp54K2XJN9MBqBJHz0EMPUVpayrRp0xg5cmTIltu+fXseeOABAO69915f256JO08//TR79uyhT58+3HXXXSFbbnJyMr/73e8AeOSRR/jiiy9CtmwT5WrKDtH+iuYawY4dOzQhIUGTkpKC6iCuzcmTJ7Vdu3YK6Lx580K+fBPdSkpKtG3btgro3LlzQ778qqoqHT58uF1b0ERhNYLImDVrFlVVVcyYMYMePXqEfPkXXXQR999/PwC/+c1vrFYQZ5599lmOHDlCbm4uU6ZMCfnyRYTf/OY3ADzxxBOcOnUq5Osw0ccSQQh9+umnvPTSSyQkJPh31uHw3e9+l8zMTNatW8fixYvDth4TXc6ePctvf/tbAH75y18iImFZz9VXX83w4cMpLi7mT3/6U1jWYaKLJYIQevLJJ6moqGD69Ol07949bOtJS0vjpz/9KYB/x2Cavr/97W8cOHCAfv36cf3114dtPSLCgw8+CMBTTz3lv421abosEYTI8ePHeeaZZwC45557wr6+73//+7Ro0YK8vDw2btwY9vUZb6kq//M//wP4Tu9MSAjvpjtp0iT69u1LYWEhr776aljXZbxniSBEXnjhBUpKShg7diwDBgwI+/pat27NnXfeCfhqIqZpW7JkCVu3bqV9+/bceuutYV+fiPDzn/8cgMcff9z6opo4SwQhUFlZye9//3sAfvKTn0RsvXfffTciwt/+9jcOHz4csfWayHviiScA+MEPfkBKSkpE1nnrrbfSpk0bNmzYwIoVKyKyTuMNSwQhsHDhQj755BO6devG5MmTI7be7t27c91111FWVsYLL7wQsfWayNq7dy8LFiygWbNmIb1uoC6pqal873vfA+CPf/xjxNZrIs8SQQg8/fTTANx1111hb7ut7vvf/z7gO63QHjvYND333HOoKtOmTSPSj2mdOXMmCQkJvP766xw9ejSi6zaRY4kgSIWFhSxYsIDk5GRuu+22iK9/woQJdOnShb1799qppE1QeXm5v7YXydqAq1OnTnz1q1+lvLycF198MeLrN5ERkkQgIhNFZLuI7BKR806gF5HHRWSj89ohIscDxlUGjJsfingi6c9//jNVVVV8/etfj/jRGkBiYqK/0/i5556L+PpNeC1YsIDPPvuMPn36hPR2JQ0xc+ZMAJ5//nnrNG6qarrcuCEvIBHfIyi7AinAJqDvBaa/G3gh4POphq4zWm4xUVlZqV26dFFA8/LyPItj//79CmhycrIePXrUszhM6F133XUK6G9/+1vPYqioqNBLL7005DdRNJFHGG8xMRjYpap7VLUMeAWYeoHpb8H3jOOY99577/Hpp5/SuXNnxo4d61kcWVlZTJo0ifLycv7+9ybx0xrg8OHDLFy4kKSkJL797W97FkdiYiIzZswAsOahJioUiSALKAz4vN8Zdh4R6QLkAEsDBqeKSIGI5IvI12pbiYjMdKYriJZOK3ej+Na3vhXxTuLqbENtev72t79RWVnJpEmTaNu2raexuOXr1Vdf5cyZM57GYkIv0nuv6cBsVa0MGNZFVXOBW4EnRKRbTTOq6jOqmququV60xVdXUlLC7NmzATw9WnNNnTqViy++mA8++ICtW7d6HY4JATepuzthL/Xp04errrqKL774grlz53odjgmxUCSCA0CngM8dnWE1mU61ZiFVPeD83QO8CwwMQUxhN3/+fE6dOsXQoUPp2bOn1+GQmprKzTffDMDLL7/scTQmWFu3bmXTpk2kp6dz3XXXeR0OcO6Ax8pX0xOKRLAO6CEiOSKSgm9nf97ZPyLSG0gHVgcMSxeRZs77NsAI4KMQxBR2blt8JC73ry83lldeecXO7ohxbvm68cYbadasmcfR+Nx8880kJiayaNEiioqKvA7HhFDQiUBVK4AfAouAj4FXVXWriDwiIoE3TJ8OvKJf3kP1AQpEZBOwDJilqlGfCIqLi3n77bdJSEjwH4VHg1GjRtGhQwf27t3LmjVrvA7HNJKq+hPBLbfc4nE057Rt25Zx48ZRUVHB66+/7nU4JoRC0kegqgtUtaeqdlPVR51hv1LV+QHTPKSq91ebb5Wq9lfVK5y/z4cinnB74403KC8vZ+zYsbRr187rcPwSExP9iemVV17xOBrTWAUFBezZs4f27dszZswYr8P5Ejcx2dlpTYtdWdwIr732GgDTp0/3OJLzuRvq7Nmz7ZYTMcq97fO0adNITEz0OJovu+GGG0hOTmbFihUcOnTI63BMiFgiaKDPP/+cJUuWkJSUxNe//nWvwznPVVddRXZ2NgcOHGDVqlVeh2MaSFX9ieAb3/iGx9Gc7+KLL2bSpElUVlYyZ84cr8MxIWKJoIHmzZtHZWUl48aNIz093etwziMi3HjjjQD+01tN7CgoKGDfvn20b9+eYcOGeR1OjaZNmwZY+WpKLBE0kNtJ5m4M0ciN7fXXX7fmoRjjlq8bb7zR84sUa3P99deTnJzM8uXL7Y6kTUR0lrQodfz4cRYvXkxiYiJf+1qtF0F7bvDgwXTs2JH9+/ezdu1ar8Mx9aSq/qPsaD7QaN26Nddeey1VVVXWPNREWCJogH/+85+Ul5czZswY2rRp43U4tUpISPD3X9iGGjs2b97M7t27yczM9OxOo/XlNj9a+WoaLBE0wLx58wCiujbgcmN0YzbRz/1fTZkyJerOFqru+uuvJyEhgaVLl/LFF194HY4JkiWCeiotLWXhwoWA774+0W7UqFGkp6ezfft2tm3b5nU4ph7ce/jEwoFGZmYmI0aMoKysjLffftvrcEyQLBHU05IlSygpKWHgwIF07tzZ63DqlJSU5L9HjdUKol9hYSHr16+nRYsWjBs3zutw6sU9ILLyFfssEdTT/Pm+i6Rj4WjNZc1DsePNN98EYOLEiTRv3tzjaOrHTQQLFiygvLzc42hMMCwR1ENVVZV/Q73++us9jqb+JkyYQFJSEvn5+Rw5csTrcMwFuAcasVS+unfvTu/evTl+/DgrVqzwOhwTBEsE9bBhwwYOHTpEx44dGTBggNfh1NtFF13E+PHjUVUWLFjgdTimFidPnmTZsmWICJMnT/Y6nAZxawX//Oc/PY7EBMMSQT24hfy6665DRDyOpmHcfgLbUKPX4sWLKSsrY+jQoUTDQ5cawspX02CJoB4CE0GscY8wFy1aRFlZmcfRmJrEcvkaOnQoGRkZ7Ny5kx07dngdjmkkSwR1OHToEAUFBTRv3tzTB9Q3VnZ2Nv369ePUqVMsX77c63BMNVVVVbz11ltAbCaCpKQkJk2aBJzr8DaxxxJBHdxzpMeOHRszZ3NU5+5g3OsgTPTYsGEDhw8fplOnTvTv39/rcBrFylfsC0kiEJGJIrJdRHaJyP01jL9NRI6KyEbndWfAuBkistN5ef+U7mrcTtavfvWrHkfSeO4Rm3UYRx/3fzJp0qSY639yTZgwgYSEBN577z1OnjzpdTimEYJOBCKSCPwBmAT0BW4Rkb41TPoPVR3gvJ5z5s0Afg0MAQYDvxaRqLm3c0VFBXl5ecC5nWksGjZsGK1atWL79u3s3bvX63BMAPcoOpbLV0ZGBkOGDKG8vJylS5d6HY5phFDUCAYDu1R1j6qWAa8A9b0Hw1eAPFUtVtVjQB4wMQQxhcTq1as5ceIEvXr1Iicnx+twGi05OZlrr70WsOp7NCkuLmbNmjUkJyfHzNXEtXETmZWv2BSKRJAFFAZ83u8Mq+5GEflQRGaLSKcGzouIzBSRAhEpiNQ90JvC0ZrLNtTo884771BVVcWoUaNo2bKl1+EEJbB8qarH0ZiGilRn8ZtAtqpeju+o/8WGLkBVn1HVXFXNjdS51osWLQJ8l/3HOvc7LFu2jLNnz3ocjYGmVb4GDRpEmzZt2Ldvn93kMAaFIhEcADoFfO7oDPNT1SJVdfc+zwFX1nderxw+fJj169eTmprK6NGjvQ4naFlZWfTv35+SkhJ7lnEUUNUmlQgSEhL4yle+ApxLcCZ2hCIRrAN6iEiOiKQA04H5gROISPuAj1OAj533i4AJIpLudBJPcIZ5zu0kHjNmTMyeNlrdhAkTANtQo8GWLVs4dOgQ7du3p1+/fl6HExKWCGJX0IlAVSuAH+LbgX8MvKqqW0XkERGZ4kz2IxHZKiKbgB8BtznzFgP/ji+ZrAMecYZ57p133gHO7TybAndDdb+b8Y67s5wwYULMnjZanXtCwvLlyyktLfU4GtMQEosdO7m5uVpQUBC25VdVVdGhQwcOHz7Mli1buOyyy8K2rkg6c+YMGRkZlJaW8tlnn9GuXTuvQ4pbEyZMIC8vj5dffplbb73V63BCZsCAAWzatIm8vDzGjx/vdTimGhH5QFVzqw+3K4trsHnzZg4fPkxWVhZ9+9Z0SURsat68OWPGjAHONX2ZyDtz5gzvvfcecO4ouqmw5qHYZImgBm7TybXXXttkqu0ut6nLEoF3VqxYwdmzZxk4cGDM3W20Lla+YpMlghosXrwYaHpHa3DuOy1evNjO9/aIW76aUv+Ta8SIEaSmprJp0yZ7GFIMsURQTWlpqb/aHutXe9akX79+tGvXjoMHD/Lxxx/XPYMJOfdouSm2oaempjJq1CjA95xvExssEVSzatUqSktLufzyy5tkZ6qI+HdA7pGpiZyjR4+yceNGUlNTGTlypNfhhIWVr9hjiaAa92itKTYLudwN1dpxI889Sh45ciSpqakeRxMe7raTl5dnzY8xwhJBNe6G2hSr7a7A873Ly8s9jia+NOX+J9cVV1zBJZdcQmFhITt37vQ6HFMPlggCHDt2jIKCApKTk/3tnE1RVlYWvXr14uTJk6xdu9brcOKGqvoTQVPsf3IlJCT4v5/1E8QGSwQB3n33XVSVYcOGkZaW5nU4YeVuqHb/+MjZu3cvn376Kenp6QwYMMDrcMLKEkFssUQQwC20TflozWUbauS5v/U111xDYmKix9GEl1u+li1bRlVVlcfRmLpYIgjgHh3H4kPqG+rqq69GRFi9ejWnT5/2Opy4EE8HGl27dqVz584UFxezadMmr8MxdbBE4HDPq09LS2Pw4MFehxN2GRkZDBo0iLKyMlauXOl1OE2eqvoPNOIhEYiI1TpjiCUCh7uRjh49mpSUFI+jiQzbUCNny5YtHD16lKysLHr27Ol1OBFh5St2WCJwLFu2DIiPZiGX+13d727CJ7B8NbX7V9XmmmuuAXz3VrLTlKObJQJHPPUPuEaMGEFSUhIFBQWcOHHC63CaNLd8uTvHeNChQwd69+5NSUkJ69at8zoccwEhSQQiMlFEtovILhG5v4bxPxORj5yH1y8RkS4B4ypFZKPzml993kjYu3cvn3zyCa1bt+aKK67wIgRPXHTRRQwZMoSqqipWrFjhdThNVmVlJcuXLwfi60ADziU+q3VGt6ATgYgkAn8AJgF9gVtEpPpN/DcAuc7D62cD/xUw7oyqDnBeU/CAW0ivvvrqJn9aX3XuhmrXE4TPxo0bOX78ODk5OXTp0qXuGZoQN/FZ+YpuoagRDAZ2qeoeVS0DXgGmBk6gqstU1T1HMR/fQ+qjhpsI4qna7rIjtvCL5/J19dVXA+du5miiUygSQRZQGPB5vzOsNncACwM+p4pIgYjki8jXaptJRGY60xUcPXo0uIgDqGpcb6jDhg0jJSWFTZs2UVRU5HU4TVI8l682bdrQv39/SktLyc/P9zocU4uIdhaLyDeBXOC3AYO7OM/QvBV4QkS61TSvqj6jqrmqmhvKpzrt2rWLAwcO0KZNmybzbOKGaN68OcOHD0dV/c9hMKFTUVHh73+Jx0QA5763209iok8oEsEBoFPA547OsC8RkfHAg8AUVT3rDlfVA87fPcC7wMAQxFRv7777LgBjxowhISE+T6Jyq+/ub2FCZ/369Zw8eZIePXqQlXWhinLT5ZYva36MXqHY860DeohIjoikANOBL539IyIDgafxJYEjAcPTRaSZ874NMAL4KAQx1Vs8V9tdtqGGT+CJCPFqzJgxiAj5+fnWTxClgk4EqloB/BBYBHwMvKqqW0XkERFxzwL6LXAR8Fq100T7AAUisglYBsxS1YglAlX1HwXH84Y6ZMgQUlNT2bx5M59//rnX4TQpbvmK5wONjIwMLr/8cs6ePWv9BFEqJG0hqrpAVXuqajdVfdQZ9itVne+8H6+q7aqfJqqqq1S1v6pe4fx9PhTx1NeOHTs4dOgQmZmZ9O1b/YzX+JGamsqwYcMAa8cNpfLycn//wJgxYzyOxlt2dlp0i89GcUdgbSBeLvuvjTUPhd4HH3xASUkJPXv2pEOHDl6H4ykrX9EtrhOBe/Qbz81CLvc3sBpB6Fiz4zmjR49GRFizZg1nzpzxOhxTTdwmAusf+LLBgweTmprKli1brJ8gROxA45z09HSuuOIKysrKWLNmjdfhmGriNhHs3LnT3z/Qp08fr8PxXGA/gV1PELzy8nLef/99wPoHXO7vYKcpR5+4TQTu0Zp7apux6wlCaf369Zw6dcr6BwJY+YpecZsIrFnofO4Rm/UTBC/wQMP4uP0Edj1B9InLRKCqtqHWYMiQITRr1owPP/zQ7jsUpMAr1o1PRkYG/fv3t+sJolBcJoLdu3dz4MABLrnkkri+fqC61NRUhg4dCmDPJwhCRUWFv3/AapxfZrXO6BSXiSCwNhCv9xeqjXXoBW/Dhg2cPHmSbt26xe39hWpjpylHp7jcC1qzUO3siC14Vr5qN3r0aABWr17N2bNn65jaRIolAvMlQ4cOJTk5mU2bNnHs2DGvw4lJVr5q597uvbS01J5jHEXiLhF88skn7Nu3j/T0dPr37+91OFGnRYsWDBkyBFX1t3Ob+qusrLT7C9XBap3RJ+4Sgdv2PXr0aOsfqIVtqI334YcfcuLECbKzs+Pu+cT1ZdcTRJ+42xO6Oze3rdKczxJB41mzUN3cbW/VqlWUl5d7HI2BOE4EtqHWbvjw4SQlJbF+/Xq++OILr8OJKVa+6tauXTt69erF6dOn+eCDD7wOxxBniaCwsJC9e/fSqlUrBgwY4HU4USstLY3c3FyqqqpYuXKl1+HEjKqqKv99miwRXJjVOqNLSBKBiEwUke0isktE7q9hfDMR+Yczfo2IZAeMe8AZvl1EvhKKeGrjFrqRI0eSmJgYzlXFPLf6bhtq/W3dupXi4mKysrLIycnxOpyoZokgugSdCEQkEfgDMAnoC9wiItUv170DOKaq3YHHgcecefvie8bxZcBE4P9zlhcWVm2vP9tQG85uZFh/bvl6//33qaio8DgaE4oawWBgl6ruUdUy4BVgarVppgIvOu9nA+PEt6VMBV5R1bOquhfY5SwvLDZu3AhYIqiPkSNHkpCQQEFBASUlJV6HExPsQKP+srKy6NatGydPnvRvl+bCTpw4gaqGZdmhSARZQGHA5/3OsBqncR52fwK4pJ7zAiAiM0WkQEQKjh492qhA8/PzWb9+PYMGDWrU/PGkVatWDBw4kIqKClavXu11OFFPVa1/oIGs+bFhpk2bRlZWVlhu2BczncWq+oyq5qpqbmZmZqOWkZiYyMCBA0lOTg5xdE2TNQ/V3/bt2zly5Ajt2rWjZ8+eXocTE9zyZQ9Cqlt5eTmrVq3i0KFDZGdnh3z5oUgEB4BOAZ87OsNqnEZEkoCLgaJ6zms8Yomg/qx/oOHc8rVixQqqqqo8jia6ffDBB5w+fZpevXpx6aWXhnz5oUgE64AeIpIjIin4On/nV5tmPjDDeT8NWKq+xq75wHTnrKIcoAewNgQxmRAYNWqUPXC8nqx/oOGys7Pp3Lkzx44dY/PmzV6HE9XCXb6CTgROm/8PgUXAx8CrqrpVRB4RkSnOZM8Dl4jILuBnwP3OvFuBV4GPgLeBf1XVymBjMqGRnp7O5Zdfbg8cr4M96KjxrNZZP1GfCABUdYGq9lTVbqr6qDPsV6o633lfqqo3qWp3VR2sqnsC5n3Uma+Xqi4MRTwmdKwdt267d+/m4MGDtGnTxh501ECWCOoW+KCjqE4Epulyz+ywG4TVzt2JuU1ppv4CE4H1E9Rs48aNYX/QkSUCc0H2IJG6uYnAHkvZcN26daNDhw4UFRXx0UcfeR1OVIrE868tEZgLyszMtAeJ1MH6BxpPRKx5qA6RKF+WCEydbEOtnT3oKHhWvmoXqQcdWSIwdbINtXZutX3UqFH2oKNGCixf4bqFQqxyH3TUpUuXsD7oyEquqZPbT7By5Up7kEg11iwUvF69etGuXTuOHDnC9u3bvQ4nqkSqfFkiMHW69NJL6d27N6dPn6agoMDrcKKKJYLgWT9B7SwRmKhiG+r59u3bx969e7n44ovtQUdBcs+4svJ1TuCDjq655pqwrssSgakXe+D4+QKvH7AHHQUnsHxZP4HPli1bKC4upnPnzmG50VwgSwSmXtwagfUTnOMmRbt+IHi9e/embdu2HDp0iJ07d3odTlQIvH4g3BcqWiIw9dK+fXt69erFqVOnWL9+vdfhRAVLBKET2E9gtU6fSJYvSwSm3mxDPaewsJA9e/bQqlUr6x8IEeuHOiewfyASJyJYIjD1Zv0E57g7q5EjR1r/QIi45WvZsmVx30+wdetWioqKyMrKomvXrmFfnyUCU2+BDxyP936CZcuWAeE/myOe9O3bl8zMTOsn4MvlKxI3MrREYOqtQ4cO9OzZ0/oJsP6BcBARq3U63O8fqQMNSwSmQdyC6R6xxKN9+/b5+wcGDhzodThNiiUCX/9ApO9oG1QiEJEMEckTkZ3O3/QaphkgIqtFZKuIfCgi3wgY92cR2SsiG52X9bpFOdtQz3330aNHW/9AiFk/AWzevJni4mI6depETk5ORNYZbI3gfmCJqvYAljifqzsNfFtVLwMmAk+ISOuA8feo6gDntTHIeEyYWT+BNQuFU58+fWjbti2fffYZO3bs8DocTwQ2C0XqQUfBJoKpwIvO+xeBr1WfQFV3qOpO5/1B4AiQGeR6jUfat29Pnz59KCkpidvnE7jNYpYIQi+wnyBemx+9KF/BJoJ2qnrIef8Z0O5CE4vIYCAF2B0w+FGnyehxEWl2gXlnikiBiBQcPXo0yLBNMOJ5Q927dy+ffPIJrUaXOdAAABr+SURBVFu3tusHwiSe+6EqKyv9/QORPCOtzkQgIotFZEsNr6mB06mvQa/WRj0RaQ/8BfiOqroPJ30A6A1cBWQA99U2v6o+o6q5qpqbmWkVCi+5BXTp0qUeRxJ57s5pzJgx1j8QJmPHjgXis59g48aNHD9+nJycnLDfXyhQUl0TqOr42saJyGERaa+qh5wd/ZFapmsFvAU8qKr5Act2axNnReT/gP/XoOiNJ9wawapVqzh79izNmtVakWty3ETg7qxM6PXo0YMOHTpw8OBBtm7dSr9+/bwOKWK8uj4l2Kah+cAM5/0MYF71CUQkBZgDvKSqs6uNa+/8FXz9C1uCjMdEQGZmJv3796e0tJT8/Py6Z2giVNUuJIsAEfEn2nirdbrfN9IHGsEmglnAtSKyExjvfEZEckXkOWeam4HRwG01nCb6sohsBjYDbYDfBBmPiZB4bB7auXMnBw4coE2bNlx22WVeh9OkxWM/QXl5uf/5xDFVI1DVIlUdp6o9VHW8qhY7wwtU9U7n/V9VNTngFFH/aaKqOlZV+6tqP1X9pqqeCv4rmUiIxyM297teffXV9nziMHPL17vvvktlZaXH0URGQUEBp06domfPnnTo0CGi67bSbBplzJgxJCQkkJ+fT0lJidfhRISbCMaNG+dxJE1fdnY2OTk5HD9+nA0bNngdTkQsWbIE8KZ8WSIwjdK6dWuuvPJKKioq/NXZpqyqqsoSQYS5v7O7g2zqLBGYmBRPG+qHH35IUVERnTp1onv37l6HExfiqXydPn2aVatWISKenIhgicA0WjxtqIFHa5G67D/euf0E77//PmfPnvU4mvBatWoVZWVlDBw4kIyMjIiv3xKBabQRI0aQkpLCxo0bKSoq8jqcsPKy2h6v2rZtS//+/Tlz5kyTP03Z6/JlicA0WvPmzRkxYgSq2qTPHiorK/Nf9m8XkkWWu2PMy8vzOJLwcr+fV+XLEoEJyvjxvgvPm3LzUH5+PqdPn6Zv374RP60v3sVD82NxcTHr168nJSWFUaNGeRKDJQITlGuvvRZo2kdsixcvBs59VxM5Y8aMISkpibVr13L8+HGvwwmLpUuXoqoMHz6ctLQ0T2KwRGCCMmjQIFq3bs2ePXvYs2eP1+GEhZvk3NqPiZyWLVsydOhQqqqqmuzDkNwDDS/LlyUCE5TExER/u6ZboJuSEydOsHbtWpKSkvwP5TGR1dRrne738rLGaYnABK0pb6jLli2jqqqKoUOH0rJlS6/DiUtNuXzt3buXPXv2+C/Q9IolAhM0d0NdsmRJk7svzDvvvANYs5CXrrrqKlq1asXOnTv59NNPvQ4npNzyNXbsWE+fb2GJwAStW7dudO3alWPHjlFQUOB1OCHlbqhf+cpXPI4kfiUlJfnPHnL/H01FtJQvSwQmJNyC3JQ21N27d7N7925at27NVVdd5XU4cc0tX4sWLfI4ktCpqKjwnxY7YcIET2OxRGBCwi3ITWlDDWwWssdSesstX0uWLKGiosLjaEJj7dq1nDhxgp49e0b0sZQ1CSoRiEiGiOSJyE7nb3ot01UGPJRmfsDwHBFZIyK7ROQfztPMTAxy2zjz8/M5ceKE1+GEhJvUvK62G8jJyaFHjx4cP36cdevWeR1OSERT+Qq2RnA/sERVewBLnM81ORPwUJopAcMfAx5X1e7AMeCOIOMxHmnVqhXDhw+nsrKySZxGWl5eHjXVduPj7jAXLlzocSSh8fbbbwPRUb6CTQRTgRed9y/ie+5wvTjPKR4LuM8xbtD8JvpMnDgROFfAY9nKlSs5deoUffv2pXPnzl6HYziXCJpC+Tp69Cjr1q0jJSUlKp5/HWwiaKeqh5z3nwHtapkuVUQKRCRfRNyd/SXAcVV1G/z2A1m1rUhEZjrLKDh69GiQYZtwmDRpEuA7YlNVj6MJjnvU6X4n471rrrmGlJQUCgoKiPV9wDvvvIOqMnr0aM9uKxGozkQgIotFZEsNr6mB06lvy69t6++iqrnArcATItKtoYGq6jOqmququZmZmQ2d3UTAgAEDuPTSSzlw4ABbtmzxOpygWCKIPmlpaYwZMwZVjfmTEqKtfNWZCJyH0ver4TUPOCwi7QGcv0dqWcYB5+8e4F1gIFAEtBaRJGeyjsCBoL+R8YyI+JuHYrkdd//+/WzevJm0tDRGjhzpdTgmQGCtM1ZVVVX5E1nMJII6zAdmOO9nAPOqTyAi6SLSzHnfBhgBfOTUIJYB0y40v4ktTWFDddugx40bR7NmzTyOxgRyy9eiRYti9ir2goICPv/8c7p06ULv3r29DgcIPhHMAq4VkZ3AeOczIpIrIs850/QBCkRkE74d/yxV/cgZdx/wMxHZha/P4Pkg4zEeu/baa0lMTOT999+P2dsGv/XWW0D0HK2Zc3r16kVOTg5FRUWsXbvW63AaJbB8RctjT4NKBKpapKrjVLWH04RU7AwvUNU7nferVLW/ql7h/H0+YP49qjpYVbur6k2q2rQfTBoH0tPTGTlyJBUVFTF5lXFpaan/5mbXXXedx9GY6kTE/3/55z//6XE0jePGff3113scyTl2ZbEJuVjeUJcvX05JSQkDBgygY8eOXodjahDL5evgwYOsX7+e5s2bR8Vpoy5LBCbk3A11wYIFMdeO6+5crDYQvcaMGUNaWhoffvghhYWFXofTIAsWLAB8ty1p3ry5x9GcY4nAhFyvXr3o1q0bRUVFrF692utw6k1VmT/fdwcUSwTRq1mzZv6rcd3/V6yYN893Pky0lS9LBCbkRIQpU3x3EnnzzTc9jqb+Nm/ezL59+2jXrp3dbTTKueUrlpqHTp8+zZIlS77UzxEtLBGYsHA31Dlz5sTMVcZz584FfJ14CQm2aUSzyZMnk5CQwJIlS/jiiy+8Dqde8vLyOHPmDLm5uXTo0MHrcL7ESrsJi5EjR5KRkcHOnTvZtm2b1+HUi5sIvvY1u+VVtMvMzGTEiBGUl5fHzDUr0Vy+LBGYsEhKSvJXf9120Wi2b98+NmzYQFpamv9pWCa6TZ3qu8tNLJSviooKfzOpG3c0sURgwsY98pkzZ47HkdTN3ZlMnDiR1NRUj6Mx9eHuUN966y3Kyso8jubCVq1aRVFREd27d6dv375eh3MeSwQmbCZMmEDz5s1Zu3Yt+/fv9zqcC3r99deB6Ky2m5p1796dfv368cUXX/ifHRGtAstXtFxNHMgSgQmbtLQ0/20a3njjDY+jqd3hw4dZsWIFycnJUXW1p6nbtGm+W5W5O9poVFVV5S//brzRxhKBCSu34M+ePbuOKb0zd+5cqqqqmDBhAhdffLHX4ZgGcMvXnDlzKC8v9ziamrk14k6dOjF48GCvw6mRJQITVpMnT6ZZs2a8//77HDp0qO4ZPOAmqWg9WjO169u3L71796a4uJjly5d7HU6N3PJ14403RmWzEFgiMGHWqlUrJk6ciKpGZa3g6NGjLFmyhJSUFP+1DyZ2iAg33XQTAK+88orH0ZyvqqqKV199FRGJ6gMNSwQm7NwN9R//+IfHkZxv9uzZqCrjxo0jIyPD63BMI7jla86cOVF39lB+fj6FhYVkZWUxbNgwr8OplSUCE3ZTp06lefPmrFy5kn379nkdzpf8/e9/B+CWW27xOBLTWP3796dfv34UFxf7byEeLdzyNX369Ki+Wj16IzNNxkUXXeQ/Gyeaqu+FhYWsWLGC1NRUO200xrmJ3N3xRoOKigpeffVVIPoPNIJKBCKSISJ5IrLT+ZtewzTXiMjGgFepiHzNGfdnEdkbMG5AMPGY6BWNG6rbVHXdddfRsmVLj6MxwZg+fTrgOwPs9OnTHkfjs2zZMo4cOULPnj0ZOHCg1+FcULA1gvuBJaraA1jifP4SVV2mqgNUdQAwFjgNBD666h53vKpuDDIeE6UmTZpE69at2bhxI5s3b/Y6HFSVl156CYBbb73V42hMsLp27crQoUMpKSmJmivZA8tXtJ4t5Ao2EUwFXnTevwjUVb+eBixU1ehI2SZimjVr5j9qe/HFF+uYOvzchHTJJZcwefJkr8MxITBjxgwgOsrXyZMn/ReRffvb3/Y4mroFmwjaqap7cvhnQLs6pp8OVG8beFREPhSRx0WkWW0zishMESkQkYKjR48GEbLxiruh/vWvf6WiosLTWNydxS233EJKSoqnsZjQ+MY3vkFKSgqLFy/2/JYms2fP5vTp04wePZqcnBxPY6mPOhOBiCwWkS01vL50Cz313XS+1hvPi0h7oD+wKGDwA0Bv4CogA7ivtvlV9RlVzVXV3MzMzLrCNlFoyJAh9OzZk8OHD3t66+CysjJ/td1NTib2paenM2XKFFSVP//5z57G8vzzzwOxURuAeiQCVR2vqv1qeM0DDjs7eHdHf+QCi7oZmKOq/uvAVfWQ+pwF/g+IzuuvTUiICHfccQcAL7zwgmdxzJ8/n2PHjnH55Zdz5ZVXehaHCb3bb78dgP/7v/+jqqrKkxi2b9/OypUrueiii/zXOES7YJuG5gPuIdUM4EI3Br+Fas1CAUlE8PUvbAkyHhPlvv3tb5OUlMSbb77JwYMHPYnhmWeeAeCOO+6I+k480zATJkygU6dO7Nmzh6VLl3oSw7PPPgvAzTffTKtWrTyJoaGCTQSzgGtFZCcw3vmMiOSKyHPuRCKSDXQCqt8M5GUR2QxsBtoAvwkyHhPlLr30UqZOnUplZaUntYLdu3eTl5dHamoq3/rWtyK+fhNeiYmJ3HnnnQA8/fTTEV9/aWmpv1nqrrvuivj6GyuoRKCqRao6TlV7OE1Ixc7wAlW9M2C6T1Q1S1Wrqs0/VlX7O01N31TVU8HEY2LDzJkzAd+RU6Q7jQOP1tLTz7vsxTQBt99+OwkJCcydOzfiNzp84403KCoq4oorruCqq66K6LqDYVcWm4gbP3483bt3Z9++fcyfPz9i6z19+rS/Wej73/9+xNZrIqtjx45MmTKFiooK/vSnP0VsvarKE088AfjKVyw1O1oiMBGXkJDAj370IwD/hhMJf/nLXzh27BhDhgxh6NChEVuvibyf/OQnAPzxj3+ktLQ0IuvMz89n3bp1ZGRkxFyzoyUC44nbbruNVq1asWLFCgoKCsK+vqqqKn/ScXcSpukaPXo0AwcO5OjRo7z88ssRWefjjz8O+PoGWrRoEZF1hoolAuOJli1b+jvTHnvssbCv76233mLbtm1kZ2dz4403hn19xlsiws9+9jMA/uu//ivsp5Lu3r2b1157jdTUVH7wgx+EdV3hYInAeOZHP/oRycnJvP7662zfvj1s61FVfvMb3wlpP/7xj0lOTg7bukz0uPnmm+ncuTM7duwI+/2HZs2aBfiuVO/YsWNY1xUOlgiMZzp27Mhtt92Gqvo3pHBYsmQJa9eupU2bNnz3u98N23pMdElJSeHee+8F4NFHH8V384PQKyws5MUXXyQhIYH77z/vvpsxwRKB8dR9991HYmIif/nLX9ixY0fIl6+q/OpXvwLgpz/9KWlpaSFfh4let99+O5deeikbNmxg3rwLXe/aeI8++ijl5eXcdNNN9OzZMyzrCDdLBMZT3bp14zvf+Q6VlZU8+OCDIV/+vHnzWL16NW3btuXuu+8O+fJNdGvevLm/XD3wwAMhv25l+/btPPfccyQkJPDQQw+FdNmRZInAeO6hhx6iefPmzJ49m/fffz9kyy0vL+eee+4B4Fe/+pU9fCZOzZw5k27durFt2zb/dSShcu+991JZWckdd9xB7969Q7rsSLJEYDyXlZXl32H/+Mc/prKyMiTLfeqpp9i1axeXXXaZ9Q3EsZSUFGbNmkViYiL/9m//RnFxcUiWm5eXx/z580lPT4/p2gBYIjBR4p577qFjx46sX7+eP/7xj0Evb9++ffz6178GfKcP2jMH4tuNN97IqFGjKC4u5r77ar3bfb2Vlpb6TxO977776NChQ9DL9JSqxtzryiuvVNP0zJkzRwFt2bKl7tmzp9HLqays1EmTJimg06ZNC2GEJpZ9/PHHmpycrIAuXrw4qGXde++9Cuhll12mZWVlIYow/IACrWGf6vlOvTEvSwRN17Rp0xTQoUOHNnoDe/zxxxXQjIwMPXDgQIgjNLHskUceUUDbt2+vR44cadQy8vLyFNCEhARduXJliCMML0sEJiYUFRVpx44dFdC77rqrwfO///777pPydO7cuWGI0MSy8vJyHTlypAI6fPjwBh9sfPLJJ5qWlqaAPvzww2GKMnwsEZiYsWbNGv/G9vjjj9d7vo8//lgzMjIU0AceeCCMEZpYtm/fPs3OzlZAv/Od72hVVVW95isuLtY+ffoooNdff71WVFSEOdLQs0RgYspf//pX/5H9rFmz6txYCwoKtG3btv6NtLy8PEKRmli0Zs0abd68uQI6Y8aMOmsG+/fv1/79+/v7BY4dOxahSEMrLIkAuAnYClQBuReYbiKwHdgF3B8wPAdY4wz/B5BSn/VaIogPv//971VEFNDrrrtODx8+fN40FRUV+vvf/96fNMaPH6+nTp3yIFoTa95++21t0aKFAtqnTx/9+OOPa5xu3rx5/k7mXr166b59+yIcaejUlgjEN65xRKSPkwSeBv6fqp53P2ERSQR2ANcC+4F1wC2q+pGIvAq8oaqviMifgE2qWue5g7m5uRqJWxcb77322mv8/Oc/p7CwkIsuuoixY8dyww030KJFC5YvX05eXh47d+4kISGBu+++204VNQ2ybt06br/9drZs2UJSUhITJkzgmmuuIScnhzVr1vDee++xZs0aAG644QaeffZZLrnkEo+jbjwR+UBVc88bUVN2aOgLeJdaagTAMGBRwOcHnJcAnwNJNU13oZfVCOLLnj17/KeD1vTq1KmTvvHGG16HaWJUcXGx3nnnnZqQkFBj+br44ov1ySefjMk+geqopUaQFKbEEygLKAz4vB8YAlwCHFfVioDhWbUtRERmAjMBOnfuHJ5ITVTKyclhwYIFbNu2jeXLl/Pee+9RUlLC4MGDGTBgABMmTCApKRJF2TRF6enpPPvsszz88MMsW7aMxYsX8/nnn5Odnc3kyZMZPXp0zD1opqHqbBoSkcXApTWMelBV5znTvEvtTUPTgInqPMxeRL6FLxE8BOSrandneCdgoar2qytoaxoyxpiGq61pqM7DKFUdH+S6DwCdAj53dIYVAa1FJMmpFbjDjTHGRFAk7jW0DughIjkikgJMB+Y77VXLgGnOdDOA8Nww3BhjTK2CSgQicoOI7MfX0fuWiCxyhncQkQUAztH+D4FFwMfAq6q61VnEfcDPRGQXvj6D54OJxxhjTMMFdfqoV6yPwBhjGq62PgK7DbUxxsQ5SwTGGBPnLBEYY0ycs0RgjDFxLiY7i0XkKPBpI2dvg+/WFtHG4moYi6thLK6GaapxdVHVzOoDYzIRBENECmrqNfeaxdUwFlfDWFwNE29xWdOQMcbEOUsExhgT5+IxETzjdQC1sLgaxuJqGIurYeIqrrjrIzDGGPNl8VgjMMYYE8ASgTHGxLkmlQhEZKKIbBeRXSJyfw3jm4nIP5zxa0QkO2DcA87w7SLylQjH9TMR+UhEPhSRJSLSJWBcpYhsdF7zIxzXbSJyNGD9dwaMmyEiO53XjAjH9XhATDtE5HjAuLD8XiLygogcEZEttYwXEXnKiflDERkUMC6cv1Vdcf2LE89mEVklIlcEjPvEGb5RREJ6F8d6xHW1iJwI+F/9KmDcBf//YY7rnoCYtjjlKcMZF87fq5OILHP2A1tF5Mc1TBO+MlbT8ytj8QUkAruBrkAKsAnoW22aHwB/ct5PB/7hvO/rTN8MyHGWkxjBuK4BWjjvv+/G5Xw+5eHvdRvwvzXMmwHscf6mO+/TIxVXtenvBl6IwO81GhgEbKll/FeBhfiexT0UWBPu36qecQ131wdMcuNyPn8CtPHo97oa+Gew//9Qx1Vt2uuBpRH6vdoDg5z3LYEdNWyPYStjTalGMBjYpap7VLUMeAWYWm2aqcCLzvvZwDgREWf4K6p6VlX3Aruc5UUkLlVdpqqnnY/5+J7WFm71+b1q8xUgT1WLVfUYkAdM9CiuW4C/h2jdtVLV94DiC0wyFXhJffLxPX2vPeH9reqMS1VXOeuFyJWt+vxetQmmXIY6roiULQBVPaSq6533J/E9u6X6M9zDVsaaUiLIAgoDPu/n/B/SP436HphzAt8DceozbzjjCnQHvqzvShWRAhHJF5GvhSimhsR1o1MNnS2+50o3ZN5wxoXThJYDLA0YHK7fqy61xR3O36qhqpctBd4RkQ9EZKYH8QwTkU0islBELnOGRcXvJSIt8O1MXw8YHJHfS3xN1gOBNdVGha2M1fnMYhM5IvJNIBcYEzC4i6oeEJGuwFIR2ayquyMU0pvA31X1rIjcha82NTZC666P6cBsVa0MGObl7xW1ROQafIlgZMDgkc5v1RbIE5FtzhFzJKzH9786JSJfBeYCPSK07vq4HlipqoG1h7D/XiJyEb7k8xNV/SKUy76QplQjOAB0Cvjc0RlW4zQikgRcDBTVc95wxoWIjAceBKao6ll3uKoecP7uAd7Fd6QQkbhUtSgglueAK+s7bzjjCjCdalX3MP5edakt7nD+VvUiIpfj+/9NVdUid3jAb3UEmEPomkPrpKpfqOop5/0CIFlE2hAFv5fjQmUrLL+XiCTjSwIvq+obNUwSvjIWjo4PL174ajd78DUVuJ1Ml1Wb5l/5cmfxq877y/hyZ/EeQtdZXJ+4BuLrIOtRbXg60Mx53wbYSYg6zuoZV/uA9zcA+Xquc2qvE1+68z4jUnE50/XG13knkfi9nGVmU3vn52S+3JG3Nty/VT3j6oyvz2t4teFpQMuA96uAiRGM61L3f4dvh7rP+e3q9f8PV1zO+Ivx9SOkRer3cr77S8ATF5gmbGUsZD9uNLzw9arvwLdTfdAZ9gi+o2yAVOA1Z8NYC3QNmPdBZ77twKQIx7UYOAxsdF7zneHDgc3OxrAZuCPCcf0nsNVZ/zKgd8C8tzu/4y7gO5GMy/n8EDCr2nxh+73wHR0eAsrxtcHeAXwP+J4zXoA/ODFvBnIj9FvVFddzwLGAslXgDO/q/E6bnP/xgxGO64cBZSufgERV0/8/UnE509yG7+SRwPnC/XuNxNcH8WHA/+qrkSpjdosJY4yJc02pj8AYY0wjWCIwxpg4Z4nAGGPinCUCY4yJc5YIjDEmzlkiMMaYOGeJwBhj4tz/D/YEb2cXKnZlAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "plot(results, interactive=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A faster implementation\n",
    "\n",
    "Although the code we wrote above works fine, it can be quite slow.  In particular, each call to `solve` necessitates rederiving the symbolic Jacobian, building new matrices and vectors and solver objects, using them once, and then destroying them.  To avoid this, we can create a solver object and reuse it.\n",
    "\n",
    "This is what the `solve` call does internally, only it then immediately discards all of this work.\n",
    "\n",
    "We start by creating a `NonlinearVariationalProblem` which gathers the information about the problem.  The residual, the solution variable, any boundary conditions, and so forth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "problem = NonlinearVariationalProblem(F, u_n1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we create a `NonlinearVariationalSolver`.  Here we provide the problem to be solved, and any options to the solver.\n",
    "\n",
    "In this case, we will modify the solver options used, noting that in one dimension, an LU factorisation produces no fill and is, obviously, an exact solve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "solver = NonlinearVariationalSolver(problem, solver_parameters={\"ksp_type\": \"preonly\",\n",
    "                                                                \"pc_type\": \"lu\"})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we just write the time loop as before, but instead of writing `solve(F == 0, u_n1)`, we just call the `solve` method on our `solver` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = 0\n",
    "t_end = 0.5\n",
    "while t <= t_end:\n",
    "    solver.solve()\n",
    "    u_n.assign(u_n1)\n",
    "    t += dt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1\n",
    "\n",
    "Compare the speed of the two implementation choices on a mesh with 1000 elements.\n",
    "\n",
    "- Hint: You can use the \"notebook magic\" `%%timeit` to time the execution of a notebook cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2\n",
    "\n",
    "Implement Crank-Nicolson timestepping instead of backward Euler.\n",
    "\n",
    "- Hint 1: The Crank-Nicolson scheme writes:\n",
    "\n",
    "   $$\\frac{\\partial u}{\\partial t} + G(u) = 0$$\n",
    "\n",
    "  as\n",
    "\n",
    "  $$ \\frac{u^{n+1} - u^n}{\\Delta t} + \\frac{1}{2}\\left[G(u^{n+1}) + G(u^n)\\right] = 0$$\n",
    "\n",
    "\n",
    "- Hint 2: It might be convenient to write a python function that returns $G(u)$ given a $u$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {
     "483f8c5b55a94cb5b9b3c6223be10e49": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "VBoxModel",
      "state": {
       "_dom_classes": [
        "widget-interact"
       ],
       "_model_module": "@jupyter-widgets/controls",
       "_model_module_version": "1.5.0",
       "_model_name": "VBoxModel",
       "_view_count": null,
       "_view_module": "@jupyter-widgets/controls",
       "_view_module_version": "1.5.0",
       "_view_name": "VBoxView",
       "box_style": "",
       "children": [
        "IPY_MODEL_8f43ada5d62046f7a6a50fbcc6587acf",
        "IPY_MODEL_ff45eb9314f84c0b99d6cbc4dafaeea7"
       ],
       "layout": "IPY_MODEL_c987f4ac87ac4cae9b51332042a41dd2"
      }
     },
     "8f43ada5d62046f7a6a50fbcc6587acf": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "IntSliderModel",
      "state": {
       "_dom_classes": [],
       "_model_module": "@jupyter-widgets/controls",
       "_model_module_version": "1.5.0",
       "_model_name": "IntSliderModel",
       "_view_count": null,
       "_view_module": "@jupyter-widgets/controls",
       "_view_module_version": "1.5.0",
       "_view_name": "IntSliderView",
       "continuous_update": true,
       "description": "index",
       "description_tooltip": null,
       "disabled": false,
       "layout": "IPY_MODEL_e34369af8020497289732bb687d2e8e2",
       "max": 50,
       "min": 0,
       "orientation": "horizontal",
       "readout": true,
       "readout_format": "d",
       "step": 1,
       "style": "IPY_MODEL_eb5e810f3cf7471c8f32d7e6ef99b66b",
       "value": 0
      }
     },
     "c987f4ac87ac4cae9b51332042a41dd2": {
      "model_module": "@jupyter-widgets/base",
      "model_module_version": "1.2.0",
      "model_name": "LayoutModel",
      "state": {
       "_model_module": "@jupyter-widgets/base",
       "_model_module_version": "1.2.0",
       "_model_name": "LayoutModel",
       "_view_count": null,
       "_view_module": "@jupyter-widgets/base",
       "_view_module_version": "1.2.0",
       "_view_name": "LayoutView",
       "align_content": null,
       "align_items": null,
       "align_self": null,
       "border": null,
       "bottom": null,
       "display": null,
       "flex": null,
       "flex_flow": null,
       "grid_area": null,
       "grid_auto_columns": null,
       "grid_auto_flow": null,
       "grid_auto_rows": null,
       "grid_column": null,
       "grid_gap": null,
       "grid_row": null,
       "grid_template_areas": null,
       "grid_template_columns": null,
       "grid_template_rows": null,
       "height": null,
       "justify_content": null,
       "justify_items": null,
       "left": null,
       "margin": null,
       "max_height": null,
       "max_width": null,
       "min_height": null,
       "min_width": null,
       "object_fit": null,
       "object_position": null,
       "order": null,
       "overflow": null,
       "overflow_x": null,
       "overflow_y": null,
       "padding": null,
       "right": null,
       "top": null,
       "visibility": null,
       "width": null
      }
     },
     "e34369af8020497289732bb687d2e8e2": {
      "model_module": "@jupyter-widgets/base",
      "model_module_version": "1.2.0",
      "model_name": "LayoutModel",
      "state": {
       "_model_module": "@jupyter-widgets/base",
       "_model_module_version": "1.2.0",
       "_model_name": "LayoutModel",
       "_view_count": null,
       "_view_module": "@jupyter-widgets/base",
       "_view_module_version": "1.2.0",
       "_view_name": "LayoutView",
       "align_content": null,
       "align_items": null,
       "align_self": null,
       "border": null,
       "bottom": null,
       "display": null,
       "flex": null,
       "flex_flow": null,
       "grid_area": null,
       "grid_auto_columns": null,
       "grid_auto_flow": null,
       "grid_auto_rows": null,
       "grid_column": null,
       "grid_gap": null,
       "grid_row": null,
       "grid_template_areas": null,
       "grid_template_columns": null,
       "grid_template_rows": null,
       "height": null,
       "justify_content": null,
       "justify_items": null,
       "left": null,
       "margin": null,
       "max_height": null,
       "max_width": null,
       "min_height": null,
       "min_width": null,
       "object_fit": null,
       "object_position": null,
       "order": null,
       "overflow": null,
       "overflow_x": null,
       "overflow_y": null,
       "padding": null,
       "right": null,
       "top": null,
       "visibility": null,
       "width": null
      }
     },
     "eaddcef7c5f24f0bae9a44e7536afc39": {
      "model_module": "@jupyter-widgets/base",
      "model_module_version": "1.2.0",
      "model_name": "LayoutModel",
      "state": {
       "_model_module": "@jupyter-widgets/base",
       "_model_module_version": "1.2.0",
       "_model_name": "LayoutModel",
       "_view_count": null,
       "_view_module": "@jupyter-widgets/base",
       "_view_module_version": "1.2.0",
       "_view_name": "LayoutView",
       "align_content": null,
       "align_items": null,
       "align_self": null,
       "border": null,
       "bottom": null,
       "display": null,
       "flex": null,
       "flex_flow": null,
       "grid_area": null,
       "grid_auto_columns": null,
       "grid_auto_flow": null,
       "grid_auto_rows": null,
       "grid_column": null,
       "grid_gap": null,
       "grid_row": null,
       "grid_template_areas": null,
       "grid_template_columns": null,
       "grid_template_rows": null,
       "height": null,
       "justify_content": null,
       "justify_items": null,
       "left": null,
       "margin": null,
       "max_height": null,
       "max_width": null,
       "min_height": null,
       "min_width": null,
       "object_fit": null,
       "object_position": null,
       "order": null,
       "overflow": null,
       "overflow_x": null,
       "overflow_y": null,
       "padding": null,
       "right": null,
       "top": null,
       "visibility": null,
       "width": null
      }
     },
     "eb5e810f3cf7471c8f32d7e6ef99b66b": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "SliderStyleModel",
      "state": {
       "_model_module": "@jupyter-widgets/controls",
       "_model_module_version": "1.5.0",
       "_model_name": "SliderStyleModel",
       "_view_count": null,
       "_view_module": "@jupyter-widgets/base",
       "_view_module_version": "1.2.0",
       "_view_name": "StyleView",
       "description_width": "",
       "handle_color": null
      }
     },
     "ff45eb9314f84c0b99d6cbc4dafaeea7": {
      "model_module": "@jupyter-widgets/output",
      "model_module_version": "1.0.0",
      "model_name": "OutputModel",
      "state": {
       "_dom_classes": [],
       "_model_module": "@jupyter-widgets/output",
       "_model_module_version": "1.0.0",
       "_model_name": "OutputModel",
       "_view_count": null,
       "_view_module": "@jupyter-widgets/output",
       "_view_module_version": "1.0.0",
       "_view_name": "OutputView",
       "layout": "IPY_MODEL_eaddcef7c5f24f0bae9a44e7536afc39",
       "msg_id": "",
       "outputs": []
      }
     }
    },
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
